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FOREWORD 


This research work deals with the normal and radial impact of composites 
with embedded penny-shaped cracks which represents a portion of the program sup- 
ported by the NASA-Lewis Research Center in Cleveland, Ohio. The program covers 
the period from February 13, 1978 to February 12, 1979 under Grant NSG 3179 and 
is conducted by the Institute of Fracture and Solid Mechanics at Lehigh University. 

Professor George C. Si h served as the Principal Investigator while Dr. E. P. 
Chen was the Associate Investigator who is now employed by the Sandia Laboratory 
in New Mexico. The capable guidance of Dr. Christos C. Chamis who acted as the 
NASA Project Manager is very much appreciated. His encouragement has led to the 
success of this work. 
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NORMAL AND RADIAL IMPACT OF COMPOSITES 
WITH EMBEDDED PENNY-SHAPED CRACKS 


by 


G. C. Sih 

Institute of Fracture and Solid Mechanics 
Lehigh University 
Bethlehem, Pennsylvania 18015 

and 

E. P. Chen 
Sandia Laboratories 
Albuquerque, New Mexico 87115 


ABSTRACT 


A method is developed for the dynamic stress analysis of a layered composite 
containing an embedded penny-shaped crack and subjected to normal and radial im- 
pact. The material properties of the layers are chosen such that the crack lies 
in a layer of matrix material while the surrounding material possesses the aver- 
age elastic properties of a two-phase medium consisting of a large number of 
fibers embedded in the matrix. Quantitatively, the time-dependent stresses near 
the crack border can be described by the dynamic stress intensity factors. Their 
magnitude depends on time, on the material properties of the composite and on 
the relative size of the crack compared to the composite local geometry. Re- 
sults obtained show that, for the same material properties and geometry of the 
composite, the dynamic stress intensity factors for an embedded (penny-shaped) 
crack reach their peak values within a shorter period of time and with a lower 
magnitude than the corresponding dynamic stress factors for a through-crack. 


* 

This work was completed when Dr. Chen was a faculty member at Lehigh University. 

- 1 - 



INTRODUCTION 


Advanced composite materials are multi-phased nonhomogeneous materials with 
anisotropic properties. This complicates the stress analysis for fracture, par- 
ticularly if the loading is time-dependent and the geometry involves sharp edges 
such as a crack. As a result, conventional and mathematical techniques for dy- 
namic fracture generally fail to yield accurate results. 

An effective approach for finding dynamic stresses in a nonhomogeneous com- 
posite containing a through crack has been developed [1] by utilizing both the 
Laplace and Fourier transforms. The transient boundary, symmetry and continuity 
conditions were formulated by integral representations in terms of the rectangu- 
lar Cartesian coordinates x and y and the results for the stress intensity fac- 
tors are determined numerically by solving a standard integral equation in the 
Laplace transform plane. The crack geometry was assumed to be extended infinitely 
in the z-direction or through the side wall of the composite specimen. Many of 
the failures in composites, however, were observed [2] to initiate from embedded 
mechanical imperfections such as air bubbles, voids or cavities. Hence, a more 
realistic modeling of the actual flaw geometry would be an embedded crack that 
has finite dimensions in all directions. This immediately suggests a three-di- 
mensional elastodynamic crack problem which cannot be solved effectively by ana- 
lytical means unless symmetry prevails. One approach for obtaining a solution is 
to extend the integral transform formulation for a through crack in rectangular 
coordinates [1] to that of an embedded crack in cylindrical polar coordinates. 

This necessitates the use of Hankel transforms instead of Fourier transforms. 

Although no attempt will be made to analyze the failure of the composite due 
to impact, the dynamic stress intensity factors k-j (t) and k£(t) can be readily 
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used in a given fracture criterion, say the strain energy density theory [3], 
for determining the allowable level of impact load. The new results can also 
assist the construction of composite materials for establishing impact tolerance. 
In this case, failure is assumed to initiate from a damage zone of material in 
the composite that can be approximated by an embedded crack. The time-dependent 
characteristics of the stresses for the through and embedded crack geometries are 
compared and studied for different elastic properties and dimensions of the com- 
posite. In particular, the phenomenon of elastic waves reflecting from the crack 
to the interfaces within the composite can be exhibited numerically when their 
neighboring boundaries are sufficiently close to one another. As time becomes 
very large, all of the results in this report reduce to the corresponding static 
solutions [4]. 


AXIAL SYMMETRIC DEFORMATION: PENNY-SHAPED CRACK 

Consider a penny-shaped crack of radius a that lies in a layer of material 
of thickness 2b with material properties y-j , , p^. This layer is bonded be- 

tween two media with properties yg, v 2» p 2 as 1 ’ll ustra ted in Figure 1. With 
reference to the system of coordinates (x,y,z), the z-axis coincides with the 
center of the crack and is normal to the crack situated in the xy-plane. The 
outer boundaries of the composite are assumed to be sufficiently far away from 
the crack such that the reflected waves will have a negligible influence on the 
local stresses. Only those impact loads that produce an axi symmetric wave pat- 
tern will be considered. 

For an axially symmetric deformation field, material elements are displaced 
only in the radial and axial direction and remain unchanged in the 0-direction. 
With reference to the cylindrical polar coordinates (r,0,z) in Figure 1, the 
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two nonzero displacement components can be expressed in terms of the wave poten- 


tials 4>-{r, z,t) and ^.(r,z,t) as follows: 

J J 



3r ” W 



+ 3i| j !i 

3z 3r r 


(D 


where j = 1 refers to the layer with the crack and j = 2 to the surrounding ma- 
terial. The four nontrivial stress components are given by 


3<j>. 3 iJk 

(«J - h fe 1 - W4 + M 2 +< 


r' . 3r '3r 3z 

J 


J T J 


i 34>.- ^ 

(a 0 ) . = 2y j F ( 3^ • 35^ + X f Z *3 

\J 


» 3<J> . 3i}>. $■; 

• + -r-i + -i) + X-V 2 d>. 

*1 3-7 '3-7 3y» v* / -i t 


(2) 


<*z>. 2y j 3z K 3z ' 3r 


J J 


~ 34). 3i^j ^ 84 , -j 4*4 

ri_ (2 — 1 + — ( — i + -1) 

L 3-7 v 4 - 3-7 / 3v< v 3v. r v> j 


(x rz ) - Pj C 3z (2 3r J - 3z J ) + 3r ( 3 p !L + r ^ 

J 


in which X. and p. are the Lamd constants and V 2 represents the operator 

«J \J 


v a .14 + 11- + 

3r^ r 3r 3z 2 


The governing equations can thus be obtained from the equations of motion which 
y i el d 
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( 3 ) 


3 2 d> . , 3d>. 3 2 d>. , 3 2 4. 

3r 2 r 3r dz* c^T at 2 


3 2 ^i i 3^ 3 2 ^,- i 3 2 ^< 

dr 1 r 3r r* dzr dy 3t z 


with c^. and c 2j - being the dilatational and shear wave speeds: 

A. + 2 V . ^ ^ 

'U-tV 1 ’ C2J= ^ 


(4) 


If the composite body is initially at rest, the Laplace transform of equations 
(3) further give 



Here, p is the transform variable in the Laplace transform pair: 


f*(p) = / f(t) exp(-pt)dt 
o 

( 6 ) 

?(t) = i / f*Cp) exp(pt)dp 
cm Bf 


The abbreviation Br stands for the Bromwich path of integration. Moreover, since 
the composite geometry is symmetrical about the xy-plane, it suffices to consider 
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only the solution in the upper half-space, z >_ 0. For the penny-shape crack ge- 
ometry, the Hankel transform pair [5] may be used: 


(7) 


f n (s) = J xf(x) J (sx)dx 
o n 


f(x) = J sf h Cs) J (sx)ds 
0 


where J n is the nth order Bessel function of the first kind. Applying equations 


(7) to (5), the following results are obtained: 


<J>?(r,z,p) = / [A^^(s,p)e Y]1 + A^ (s,p)e Yl1 ] J Q (rs)ds 
o 

^*(r,z,p) = / [B^Cs,p)e Yzl + B^Cs,p)e Yzi ] J^rsjds 


(8) 


for the cracked layer and 


4»2(r,z,p) = / C^(s,p)e Yl2 J Q (rs)ds 
o 


i|£(r,z,p) = / C^(.s,p)e J-j (rs)ds 


-Y 22 z 


C9) 


for the surrounding material. The quantities y. . are given by 

* J 


, , n 2 1/2 „2 1/2 

*1J - Is* ♦ fsj) . Y 2 j - (s* + frj-) 


( 10 ) 
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The six unknowns A^, A^,..., are determined from a given set of tran- 
sient boundary, symmetry and continuity conditions. 


NORMAL IMPACT 

* 

Let the penny-shaped crack be subjected to a uniform impact load such that 
the upper and lower surface will move in the opposite direction. The magnitude 
of this normal load is a Q and since it is applied suddenly from t = 0 and main- 
tained at a constant value thereafter, the Heaviside unit step function, H(t), 
will be used, i.e., -cr Q H(t). Making use of equations ( 6 ), the conditions on the 
plane z = 0 for r < a and r >. a take the forms 


(O tr,o,p) - - r 2 -; 
2 1 P 



) (r>o,p) = 0, 0<r<a 
1 


(u*) (r,o,p) = 0; (x*_) (r,o,p) = 0, r > a 
2 1 rz 1 


01 ) 


If the interfaces at z = ±b is bonded perfectly, the stresses and displacements 
can then be considered continuous across these planes, i.e.. 


(a*) (r,b,p) = (a*) (r,b,p) 

2 1 2 2 

(12) 


•ff if 

(x r _) (r,b,p) = (t ) (r,b,p) 
rz i rz 2 


if 

There is no loss in generality in formulating the problem in terms of a uniform 
step load. The principle of superposition may be used to obtain the solution for 
general loading from a series of step loading solutions as discussed in [1]. 
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and 


(uj (r,b,p) = (u ) (r,b,p) 
r 1 r 2 

03 ) 

(O (r,b,p) = (u ) (r,b,p) 
z 1 2 2 

Under these considerations, the six functions A^, A^,..., may be ex- 
pressed in terms of a single unknown A(s,p) as indicated by equations (A.l) in 
the Appendix. 

FtedhoZm lrvto.QKal e.quationi>. Without going into details, the function A(s,p) 
can be obtained from the system of dual integral equations 


/ A(s,p) J Q (rs)ds = 0, r > a 


/ sPj(s,p) A(s.p) J 0 Crs)ds = - z^dg - K.p p. r < 
in which P j (s ,p) is a known function: 


04 ) 


* j ( s » p ) . 4 U^,)‘ - s^ 21 ][* (2 > - 

+ s(s^I 1 )e" tYll+Y2,)b [y 21 (^>6< 4 > - 5 < 2 >6< 3 >) - Y„] 

+ [J (s 2 + y 2 ] ) 2 + s a Y 11 Y 21 ][« C 4 ) e’ 2Y2lb - « ( 1 ) e' 2Yllb ]) ( 15 ) 
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The form of A(s,p) that satisfies equations (14) can be found from Copson [6]: 


A(s,p) - - /J s 2n 1 pCl-xf ) { J 1/2 ( sa ^) d ? 


a o a 


5/2 


(16) 


* 

Here, Jy 2 is the half order Bessel function of the first kind and Aj(5,p) satis- 
fies the Fredholm integral equation 


At(5.p) + / A T (n,p) M T (5,n,p)dn = C (17) 

i 0 

whose kernel 

00 

Mj(.5,ri,p) = J s[Pj(— ,p) - 1] J-j^CsS) J-j ^ 2 Csn)ds 

= | 7 DMf.p) - 1] sin(s5) sin(sn)ds 08) 

o 

is symmetric in 5 and n. Figures 2 to 4 show the numerical results of equation 

(17) by varying y 2 /y.| and a/b while p-j = p 2 and v-| = v 2 = 0.29 are kept the same 

* 

for all cases. The function Aj(5,p) evaluated at the crack border, 5 = 1, gov- 

* 

erns the contribution of the geometric and material parameters on k-j(p) which 
represents the Laplace transform of the stress intensity factor. 

S&iz&A fiacAon. fan. nonmal AmpacA. In order to evaluate k-j(p) or 

k-j Ct ) , the stresses in the matrix layer are first expanded in terms of the local 
coordinates r-j and 0-j for small values of r-j . The local coordinates (r^,0-j) are 
related to (r,0) in Figure 1 as follows: 
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a + r-j cos8-| = r cose 


09 ) 


r-| sinG-j = r sine 


The leading term in the Laplace transform of the local stresses that possess the 
1M7 singularity is 


* A, (1 ,p) « 

k,(p) = — b vj a 

1 vr ' P IT 0 


(20) 


Application of the Laplace inversion theorem yields the dynamic stress field 
around the crack border as a function of time. The result is 


k,(t) e, 0 , 30, 

(a r ) ] (r 1 ,e 1 ,t) = cos ^Ml - sin j- sin + 0(rj) 

k i (t) 0-j 

= — 2v-| cos + 0(rp 


(^q) (r, >9, ,t) 

0 1 1 1 ^ 


30, 


k, (t) 0, 0, ju i 

[a) (r, ,0, ,t) = cos -J- (1 + sin sin -$*-) + 0(r?) 

z 1 1 1 t \ 


k,(t) 0 


30, 


M Vt/ Vi w>u-» 

C T rz ) = cos sin ^ cos + 0(rj) 


fir. “ 


( 21 ) 


and k-j(t) becomes 


„ 2a o^ 1 1 p pt r 

Br 


¥‘> -H-srL 


( 22 ) 


Note that equation (20) is, in fact, the Laplace transform of equation (22). 
Hence, the functional dependence of r-| and 0-j is not affected by the Laplace 
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transformation and can be evaluated separately. This observation was first made 
by Sih, Ravera and Embley [7]. 

•ft 

Making use of the results for Aj Cl »P) in Figures 2 to 4, k^(t) in equation 
(22) can be found as given in Figures 5 to 7. The dynamic stress intensity 
factors k-j (t) for the penny-shaped crack exhibit an oscillatory behavior rising 
quickly to a peak. As time increases, all curves approach the static value of 
k-j = 2cr Q >/a/'iT [4]. For a crack diameter to layer thickness ratio of a/b = 1, the 
peaks of the k^(t) curve are sensitive to changes in the shear moduli ratio y 2 /y.j . 
Figure 5 indicates that k-j (t) tends to decrease in amplitude as y 2 /y.| re duced 
from 0.1 to 10.0. The influence of the composite interface on k-j(t) is exhibited 
in Figures 6 to 7. When the shear modulus of the surrounding material y 2 is 
much smaller than the matrix layer with y.j, the dynamic crack border stress in- 
tensity increases as the crack diameter becomes large in comparison with the layer 
thickness. This effect is clearly evidenced in Figure 6. As expected, k-j(t) in- 
creases with decreasing a/b when the shear modulus of the cracked layer is made 
smaller than the surrounding material, i.e., y. 1 < u 2 as illustrated in Figure 7. 
The result of Embley and Sih [8] is recovered for the homogeneous case, y-j = y 2 « 

RADIAL IMPACT 

If the penny-shaped crack is sheared uniformly in the radial direction such 

•ft 

that axial symmetry is preserved, then <}>.(r,z,p) and i|/.(r,z,p) in equations (8) 

J J 

and (.9) remain valid. Let this shear of magnitude t q be applied suddenly and 
hence the surface tractions, -x Q H(t), are to be specified for 0£r<a with H(t) 
being the Heaviside unit step function. Laplace transform of the conditions on 
the plane z = 0 thus become 
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( 23 ) 


(t* z ) (r,o,p) = - (a*) (r,o,p) = 0, 0<r<a 


•fa ^ 

(uj (r,o,p) = 0; (ct_) (r,o,p) = 0, r > a 
r 1 z 1 

Continuity of the stresses across the interface z = b is satisfied if 


(a ) (r,b,p) = (a) (r,b,p) 
Z 1 z 2 


(<*!,) (r,b,p) = (a*) (r,b,p) 
rz 1 rz 2 


and the same requirement is imposed on the displacements: 


(24) 


(uj (r,b,p) = (u ) (r,b,p) 
r 1 r 2 

(25) 

(u ) (r,b,p) = (u ) (r,b,p) 
z 1 z 2 

X ntzgtuil equationi. As in the case of normal impact, the six unknown func- 
tions A^(s,p), A^(s,p),..., C^(s,p) in equations (8) and (9) can be ex- 
pressed in terms of a single unknown B(s,p). Refer to equations (A. 5) in the 
Appendix. Hence, equations (24) and (25) are satisfied. The remaining boundary 
conditions in equations (23) are employed to obtain the system of dual integral 
equations 
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/ B(s,p) J.|(rs)ds = o, r >_ a 


/ sP n (s,p) B(Sjp) J, (rs)ds = - ^ r < 


( 26 ) 


in which 


Pll(s.p) “ Pj (s ,p) 


(27) 


where Pj(s,p) is already known through equation (15) while Aj(s,p) and Ajj(s,p) 
are given by equations (A. 2) and (A. 6 ), respectively. 


Solving for B(s,p) [ 6 ], it can be shown that 


T a 5/2 -j 

A n«-p> 


(28) 


'K 

and Aj j (?,p) satisfies the Fredholm integral equation of the second kind: 


Ajj(?»p) + J Ajj(n.p) M n (5,n,p)dn = 5 


(29) 


whose kernel takes the form 


Mn(5,n,p) = / s C P j j (■gf» p) " J 3/2 


(30) 


Plots of Ajj(l,p) as a function of c 2 -|/pa are shown in Figures 8 to 10 for dif- 

* 

ferent values of y 2 /y.| and a/h. The curves show that Ajj(l»p) rises rapidly at 
first and then levels off. 
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Stfic&i Intensity iacton. fan. MxUaZ impact. The dynamic crack border stress 
field corresponding to radial shear can be obtained in the same way and expressed 
in terms of the coordinates (r-j , 0^ ) in equations (19): 


(&y.) (r-| »0-| >t) = 

r 1 11 


kp(t) 0, 0, 30, 

= sin (2 + cos y~ cos -jp-) + 0(rp 


ko(t) 0-t 

CO (r, ,0, ,t) = 2v, sin 2 + 0(r?) 

0 1 1 1 S2r^ 1 L 1 


k 2 (t) 0, 0, 30, 

(O (r, ,0, ,t) sin cos -rr- cos + 0(r?) 

z 1 1 1 c c L \ 


(31) 


( T rz^ r r 9 l 5t ) = cos T ^ " sin T sin + °^ r l^ 


Note that k 2 (t) can be evaluated from 


k 2 (t) = 


V® , A IlO’P) 
w J r P 


e pt dp 


(32) 


once Ajj(l,p) as given by Figures 8 to 10 is known. 

The numerical results in Figures 11 to 13 for k 2 (t) as a function of time 
refer to = p 2 and v-j = v 2 = 0.29. The curve with y-j = y 2 is the solution for 
the homogeneous material treated previously by Embley and Sih [8]. In general, 
k 2 (t) oscillates with time and can be greater or smaller than the corresponding 
homogeneous solution depending on whether y 2 /y-| < 1 or y 2 /y-j > 1. Figure 11 
displays the variations of k 2 (.t) for different values of y 2 /y-j while a/b is 
fixed at unity. The influence of the ratio of crack size with layer thickness 
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is exhibited in Figures 12 and 13 for p 2 /y.j = 0*1 and = ^.0, respectively. 

These two cases show the opposite effect which is to be expected. 

CONCLUDING REMARKS 


The previous discussion has shown that the dynamic stress intensity factors 
for an embedded crack can be evaluated analytically by a method similar to that 
developed for a through crack [1]. An important consideration is to compare the 
results for these two crack configurations and to draw some general conclusions. 
First of all, the k-j (t) or k 2 (t) factor for the penny-shaped crack tends to rise 
more quickly than the through crack, i.e., the peak value of k-j(t) or k 2 (t) is 
reached within a shorter period of time. This is because waves emanating from 
the neighboring points on the periphery of the penny-shaped crack interfere with 
each other much earlier as compared to a line (or plane) crack where the waves 
must travel from one end to the other before interference can take place. In 
general, the maximum value of k-j(t) or k 2 (.t) for an embedded crack is lower than 
that for a through crack. For example. Figure 5 gives a peak value of approxi- 
mately 1.6 for TTk-j (t)/2cr 0 »^a which corresponds to a/b = 1.0 and yg/y-j = 0*1* This 
occurs at c^t/a * 1.6 and yields k^Ct) = 1.02 o 0 ^a. The corresponding case of a 
through crack [1] renders k-j(t) = 2.40 a Q Va and c^t/a = 3.0. The difference in 
k.j(t) is more than a factor of two and is more pronounced as the ratio a/b is in- 
creased. For embedded cracks that are non-circular in shape, approximate esti- 
mates of k^(t) can be made by taking the solution for the through crack as an 
upper limit and that of the circular crack as a lower limit.- 

In the absence of axisymmetry, the dynamic stress analysis will become ex- 
ceedingly difficult and it will be more feasible to solve the crack problem nu- 
merically. In such cases, the solutions obtained here can be used to guide the 
development of numerical procedures. 
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APPENDIX: EXPRESSIONS FOR A^(s,p),..., C^(s,p) 

Normal impact. The functions A^(s,p), A^(s,p),..., C^(s,p) for the 
wave potentials in equations (8) and (9) can be expressed in terms of a single 
unknown A(s,p) for normal impact 

A^d.p) - [i (s^1,)(6 C2) ♦ - SYll e' Wll+Y2llb ] 

A< 2 >(s,p) - - [sy„e' lYll+Y2l)b <^> 


A(,s,p) 


l< 1) (s,p) = -U< 1 )A( 1 )e' Yllb + S ' 2 >A< 2 >e Yllb ] 


|( 2 >(s,p) = - [6< 3 V 1 V Yllb + , 5 ( 4 >A< 2 >e Yl l b ] 


(A.l) 


C (1) (s,p) 


Y )2 b 


[(s 2 -r 11 Y 22 )A tl) e' Yllb ♦ (s^ 11 Y 22 )A< 2 >e Yl1 ' 


s " Y 12 Y 22 

s(T 21 -T 22 )B (,, .' Y2lb ♦ s(Y 2]+ Y 22 )B( 2 )e Y2lb ] 


C (2) (s,p) -p§ 


Y 22 b 


y 12 y 22 [s(Y i2' Y ii )A(1)e Y " + ^n^’ 1 

♦ (s^-Y 21 Y 12 )B< 1 )e' Y2lb ♦ (s^ 21 Y 12 )B' 2 >e Y2lb ] 
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in which Aj stands for 


a n ^ P 2 r^C2) . ^C3) " 2 ^ Y n +Y 2i^ b 

A I (S ’ P) = 2fer Y n C5 + 5 e 


+ <5^e 


■2y 21 b 


'21 


+ «Ci>.-*n b ] 


(A. 2) 


and <5^, 6^,..., <5^ are further expressed in terms of e^\ e^,..., e^ 
as the following: 

«« 1 >(s.p) - (e^V 6 * - e< 2 V 7 >)/te ( 1 V 6) - e< 2 V 5 >) 


6< 2 >(s,p) - (e^e( 6 ) - e< 2 V 8 V(e (1 >e< 6 > - e^e^) 

(A. 3) 

S (3) (s>p) . (e O) e (7) . e C3) e (5)) /( e(l) e ( 6 ) . e (2) e (5), 

6^ (s,p) = (e^ 1 V 8 > - .We^)/(e^W - e< 2 V 5 >) 


The quantities in equations (A. 3) are complicated functions of the materials pa- 
rameters and transform variables. They are given by 

e (1) (s,p) = - sy 21 + ^ ■ (p - ^ 2 2 -- y [| (y 21 -y 22 )(s 2 +y| 2 ) + T 2 2 (s2 - y 2^12^ 

e (2) (s,p) = sy 21 - yi ( s ^Y^ 2 Y 22 ) ^ kfzi+VzzHs^h) ' T 22 (s2+Y 21 y 12^ 

e (3) (s,p) = \ (s 2 +y 21 ) - ~ (s 2 -y^ 2 y 22 ) ^ ^ s2+ ^22^ s2 " y 11 y 22^ 

+ s2 T 22 (Y ir y 12 )] 
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e (4) (s,p) = 1 (s 2+ Trl! ) - ^(s^v-g) Cs 2 +y| 2 ) C® 2 ^ ^ag ) 

” s ^22 ^11^1 2^ 

e (5) (s,p) = - 1 (s 2 +y^) + te\z^Z-\-'< 2 z' 

+ j (s 2 +y| 2 )(s 2 -y 21 y 12 )] 

(A. 4) 

e^ 6) (s,p) = - 1 Cs 2+ y|-j ) “ y 1 ts ii -Yi 2 Y 2 2) Cs2y 1 2 (y 21 +y 22 ) 

- \ (s 2 +Y22)( s2+ Y21 y 12^ 


(7) s ^2 1 

e (s,p) = s Yll - Cy^^-Yh"^ + 2 (s2+Y 22^ Y n" Y 12^ 

e^ 8 ^(s,p) = - SY n - p 1 Cs 2 -Yi2Y2 2 ) ^ Y 12^ 2+Y n Y 22^ “ \ ^ s2+Y 22^ Y n +Y 12^ 

Radial impact. For radial impact, A^(s,p), A^(s,p),..., C^(s,p) in 
equations (8) and i9) can be expressed in terms of B(s,p) as 


aW(..p) = - [ S y 21 ( 5 (2) - fi W) e -^2i b ) + 1. (s2n | l)e - (Y n +Y 2i )b ] 


(A. 5) 


A< 2 >(s, p) - [sr 21 e' 2Yllb H<’> - Zy * b ) ♦ 1 (sHy|,)e' (Yll+Y2l)b ] 


2 ■ <21 


x B(s,p) 
A II 
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where 


. _ p 2 r,(2) *C3)-“ 2 ^ Y n 4T 21^ b .(4)" 2y 21 b 

A n -257*21 [<S + 5 e " 6 e 


'21 


- 5 o) e - 2Y ” b ] 


(A.6) 


The remaining functions B^(s,p), B^(s,p), etc., can be related to B(s,p) 
through A^(s,p) and A^(s,p) since the last four expressions in equations 
(A.l) for normal impact also apply to radial impact. 
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Figure 5 - Dynamic stress intensity factor k-, ( t) for penny-shaped 
crack with a/b =1.0 
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Figure 6 - Dynamic stress intensity factor k, (t) for penny-shaped 
crack with =0.1 
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Figure 7 - Dynamic stress intensity factor k, (t) for penny-shaped 
crack with = 10*0 
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Figure 10 - Variations of Ajj(1,p) with Cg-j/pa for = 10 and varying a/b 
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Figure 11 - Stress intensity factor k 2 (t) versus time for a penny-shaped 
crack with a/b =1.0 
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Figure 13 - Stress intensity factor k~(t) versus time for a penny-shaped 
crack with yg/v-j = 10.0 
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COMPUTER PROGRAMS 


kxJjaJL impact 


3 

3 

3 

3 

3 

3 

3 

4 

5 

20 


20 

22 

23 

24 

31 

33 

34 

35 
37 
41 

46 

47 
54 

56 

57 
61 
64 
64 
72 
72 

74 

75 

102 

103 

105 

106 
111 
120 
120 
131 
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141 
144 
146 
155 
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160 
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164 
166 
202 
202 
204 
207 


* 

#• 


2 

Ki 
K2 
— K3 


PROGRAM SETA (INPUT , OUTPUT , PUNCH , PL OT ,T APE 
REAL NON (4) ♦ F (4 ,4 , i> ♦ G (4 t4> ♦ D(4> , PT (4) 
REAL B( 4) ♦ C ( 4) 

REAL LPC50) , 0TA(50) 

EQUIVALENCE (NON , B) 

CCMMON Kl» K2.K3 ♦ K4 
COMMON/ AliX/HtP t PKl,PK2,eMU,X,Y 
LP ( 1)=0 • 0 
OTA (1) = 0* 0 

READ 2, K1,K2,K3,K4 - ' - 

FORMAT ( 12) 

= ORDER CF SYSTEM OF EQUATIONS 
= NO. OF DISTINCT KERNELS 
DATA POINTS 

DATA SETS TO BE EVALUATED 
POINTS 


99 = 


= NO. OF 
K4 = NO. OF 
SET UP OATA 
AK=K3 
DO 5 N= 1, K3 
AN=N 

5 PT(N)=AN/AK 
SET UP INTEGRATION HATRIX 
M=K3-2 
N=K3-1 
A=K 3 

A=l./ (3.*A) 

DO 10 K=2,M,2 
10 0<K)=2.*A 

DO 15 K=1,N,2 
15 D(K1=4.*A 

D(K3)=A 

CALCULATE NONHOMOGENEOUS TERMS 
RHS=1.0 
DO 22 1=1, K2 
PRINT 9 
9 FORMAT (1H1) 

REAO 61, EMU 
61 FORMAT ( F10. 5) 

DO 999 11=1, K4 
DO 35 N=i,K3 
35 NON(N)=RHS*PT(N> 

CALCULATE KERNEL MATRICES 
CALL CONST ( I ) 

DO 20 N=1,K3 
00 20 M=1,K3 


25 

30 

20 


40 


999 


66 


22 


IF ( M-N) 25,30,30 
F(H,N,IJ=F (N ,M, I) 

GO TO 20 

F (M,N,I) =FU(I,PT (M),FT(N)) 
CONTINUE 

CALL CHANGE (F,G,D, I) 

CALL LINEQCG ,B, C, K3) 

DO 40 L=l, K3 
PRINT 6 ,PT( L) ,NON (L) 
FORMAT(5X,F8.4,F15.6) 

CONTINUE 

LP(II+l)=NON(K3) 

DT A (II+1)=P 
CONTINUE 

PUNCH 66, (DTA(IX) ,LP(IX) ,IX=1,19> 
FORMAT ( 2F10. 5) 

CALL LAPINV(DTA,LP) 

CONTINUE 

END 


FUNCTION SIMPtl , A ,B) 

6 COMMON/ AUX/H,P,PK1,PK2,BMU,X,Y 

6 DEL = 0. 25 4 (B-A) 

10 IF (DEL) 40,45,50 

12 45 SIMP=0. 0 

13 RETURN 

14 50 CONTINUE 

14 SA=Z(I, A)+Z(I,B) 

26 SS=Z(I, A+2.*DEL) 

35 SC=Z(I, A + OED+Z (I,A + 3.*OEL) 
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PLOT) 



53 


Sl= (OEL/3. >* (SA+2.+SB+4.*SC) 


61 


IF CS1 . EG. 0 . 0 ) GO TO 45 


62 


K=8 


63 

35 

S3=SB+SC 


65 


0£L=0.5*CEL 


67 


,SC = ZCI, A + DEL) 


75 


J=K-1 ■ - 


77 


00 5 N=3,J,2 


100 


AN=N 


101 

5 

SC=SC+ZCI, A*AN*DEL) 


113 



S2=CDEL/3.)*CSA*2.*SB+4.*SC) 


122 


0IF=ABSttS2-Sl) /SI) 


125 


ER= 0.01 


127 


IF ( DIF-ER) 30,25,25 


131 

30 

SIMP=S2 - 


133 


RETURN 


133 

25 

K=2*K 


134 


S1 = S2 


136 


IF CK-2048) 35,35,40 

w 

14 0 

40 

PRINT 42,I,A,e 


152 

42 

FORMAT t5X ,* INT. DOES NOT CONVERGE *,I3,2F9.4) 


152 


PRINT 60 , X, Y 


162 

60 

FORMATS 2F 10*5) 

• 

162 


DO 70 J=l,10 


166 


DIP=J 


167 


OIP=OIP/10. 


171 


W=Z Cl , DIF) 


175 


PRINT 60, W 


202 

70 

CONTINUE 


206 


CALL EXIT 


207 


END 



SUBROUTINE CHANGE CF, G,D, I) 


7 


REAL F(4,4,il ,GC4,4) ,0(4) 

7 


COMMON K1,K2,K3,K4 

7 


DO 10 N=1,K3 

10 


DO 10 M=i,K3 

11 


G(M,N) =F(M,N,I) *0(N) 

24 

10 

CONTINUE 

30 


DO 20 N=1,K3 

31 

20 

G(N,N)=G(N,N)*1.0 

40 


RETURN 

41 


ENO 


SUBROUTINE LINEG (A ,B,T,N) 


7 


REAL A(N,N) , E(N) ,T(N) 

7 


00 5 1=2, N 

10 

5 

A(I*1)=ACI»1)/A(i,l) 

17 


DO 10 K = 2 , N 

20 


M=K-1 

22 


00 15 1=1, N 
T C I ) =A C I , K) 

23 

15 

33 


DO 20 J = 1 ,M 

34 


At J,K) =TIJ) 

41 


J1=J+1 

43 


DO 20 I=J1, N 

44 


T 1 1 ) =T ( I) -A CI,J)* A(J,K) 

55 

20 

CONTINUE 

61 


ACK,K)=T (K) 

65 


IFCK.EQ.N) GO TO 10 

66 


M=K + 1 

70 


00 25 I=M»N 

71 

105 

25 

10 

A Cl , K)= TCI) /ACK,K) 
CONTINUE 

♦ BACK SUBSTITUTE 

110 


DO 30 1=1, N 

111 


T (I) =Bt I) 

114 


M=I + 1 

11 € 


IFIM.GT.N) GO TO 30 

121 


DO 30 J=M,N 

122 


B(J)=8< J)-ACJ,I)*T(I) 

132 

30 

CONTINUE 

136 


DO 35 1 = 1, N 
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137 

K=N + 1 -I 

141 

B(K>=T(K)/A (K,K) 

146 

Kl= K- 1 

150 

IF(Kl.EQ.O) GO TO 35 

151 

DO 35 J 1=1 * K1 

152 

J= K-Ji 

1.54 - 

T(J)=T( J)-Af J,K)*B(K) 

162 

35 CONTINUE 

167 

RETURN 

167 

END. 



6 

6 

7 

10 

11 

12 

13 

20 

21 

23 

25 

32 

33 

36 

37 
*»1 
47 
4 7 


FUNCTION FU (I, A, B) 

C0MM0N/AUX/H,P,PK1,PK2,BMU,X,Y 

X=A 

Y=B 

IF ( A*B) 5f 10 > 5 
10 FU= 0 • 0 
RETURN 

5 SUM=SIMPCI,0.0,5.Q) 

ER=0.01 . 

DEL =5. 0 - 

20 UP=DEL+5 • 0 

ADOL=SI HP ( I , DEL , UP) 

DEL =UP 

TEST=ABS(AODL/SUM> 

SUM=SUM +ADDL 
IF (TEST-ER) 15,20,20 
15 FU=SQRT (X*Y) +SUM 
RETURN 
END 


3 

3 

5 

6 
15 
24 
31 
31 

33 

34 

36 

37 

41 

42 

44 

45 
62 

62 

63 


SUBROUTINE CONST (I) 

COMMON/ AUX/H, Pi PK1,PK2, BMU, X,Y 

PR1=0.29 

PR2=0.29 

PK1 = SQRT( (l.-2.*PRl) /<2.*(1.-PR1) ) ) 

PK2=SQRT( (l.-2.*PR2) /(2.* (1.-PR2) ) ) 

READ 1,P 

1 FORMAT CF10 . 5) 

HH= 0.1 
HH=10.0 
HH=5. 0 
HH=4.0 
HH=1.0 

HH= 0. 5 
HH=2. 0 
H=1./HH 

PRINT 2,BMU,PR1,PR2,HH,F „ , 

2 FORMAT ( /////5X, * MU2/MU1 =*F6.2,* NU1 =*F4.2,* NU2 =*F4. 2///5X ,* A 
1/H =*F4.2,* C21/PA = »F*.2) 

RETURN 

END 


5 

5 

23 

25 

27 

30 

31 
31 

33 

34 
36 
40 
46 
55 
64 
73 
77 


5 

10 


FUNCTION Z(I,S) 

COM MON/ AUX/H, P,PK1,PK2»EMU,X»Y 

BESJH(A)=SQRT<2.*A/PI)*SIN(A)/A 

PI=3. 1415926 

IF(S-0. 0)5,5,10 

Z=0.0 

RETURN 

CONTINUE 


PP=P*P 
Ci=PKi*PKi 
C2=PK2*PK2 
CC— 1, — C 1 

GA=SQRT(S*S+C1/PP) 

GB=SQRT(S*S+1./PP) 

GC=SQRTfS*S+C2/BMU/PP> 

GD=SQRT (S*S+1./BMU/PP) 

AA=S*S+l./PP/2. 

A9=l.-BMt 
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100 

103 

110 

115 

123 

131 

137 

145 

152 

157 

162 

164 

171 

174 

200 

203 

206 

211 

214 

217 

222 

226 

231 

233 

235 

236 
240 
242 
244 

246 

247 
252 
256 
262 
266 
275 
301 
307 
312 
321 
332 
337 
347 
35 0 


AC=S*S-GC*GO 
AD=(GB-GD)/AC/PP/2.*BMU 
AE= (GB+GD) / AG/PP/2,*BMU 
AF = (S*S-GA*GD)/AC/PP/2.*EMU 
AG= (S*S+GA*GD) /AC/PP/2, *EMU 
AH=(S*S-GB*GC)/AC/PP/2.*BMU 
AI=(S*S+GB*GC)/AC/PP/2.*EHU 
AJ= (GA-GC) /AC/PP/ 2,* BMU 
AK= (GA+GC) /AC/PP/2.*BMU 
Ai = -(AB*GB-AD) 

A2=AB*GB-AE 

A3=AA-BMU*S*S-AF 

A4=AA-BHU*S*S-AG 

A5=-AA+ BMU*S*S+ AH 

A6=-AA+BMU*S*S+AI 

A7=S*(AE*GA-AJ> 

A8=-S*(AE*GA-AK> 

BA=A1+A6-A2*A5 

BB=A3*A6-S*A2*A7 

BC=A4*A6-S*A2*A8 

ED=S*A1*A7-A3*A5 

BE=S*A1*A8-A4*A5 

B1=BB/B A 

B2=8C/B A 
B3-BD/BA 
B4= BE/BA 
EA=2.*GA*H - 
EB=2.*G B*H 
EC=(EA+EE)/2. 

ED=2.*EC 

El-EXP (-EA) 

E2=£XP(-EB) 

E3=EXP(-EC) 

E4=EXP(-ED) 

DL=B2+B3*E4-«'B4*£2+B1*E1 

0i=2.*PP/CC/GE/DL 

D2=AA*AA-S*S*GA*GB 

D3=B2-B3*E4 

D4=2.*AA*(GB*(B1*B4-B2*B3)-S*S*GA>*E3 
D5=(AA*AA+S*S*GA*GB)» (B4*E2-SI*E1> 
F=Dl*(D2*D3«-0*H-05) 

Z= (F-S) *E£S JH(S + X) *BESJH(S*Y) 

RETURN 

ENO 


5 

5 

5 

5 

5 

16 

16 

30 

30 

34 

34 

40 

44 

44 

65 

65 

67 

73 

73 

112 

112 

116 

121 

130 

130 


SUBROUTINE LAPIN V tGL AH, PHI) 

C THIS PROGRAH EVALUATES THE COEFFICIENTS FOR SERIES 

C OF JACOBI POLYNOMIALS WHICH REPRESENTS A LAPLACE 

C INVERSION INTEGRAL 


INVERSION 
REAL MUL 
OIMENSI ON 
DIMENSION 


A (50), GLAM (50), PHI (50) ,C(4,50) 

__ 3K(101) *TT (101) 

COMMON/ 2/TI,TF,DT,MN,EK,TT 
READ 1, NN ,MN,MM 

1 FORMAT ( 312) 

READ 2, TI ,TF ,OT 

2 FORMAT (3F10,5) 

PRINT 95 

99 FORMAT ( 1H1) 

CALL SPLICE (GLAM, PHI, MM, C) 

PRINT 101 

101 F0RMAT(/////5X,* GLAM PHI *) 

PRINT 102, (GLAM (I) ,PHI(I) ,1=1, MM) 

102 FORMAT(5X,F10.5,5X,F10. 5) 

Mll=MM-i 

PRINT 300 

300 FOR MAT ( ////5X,* C(1,J) C(2,J) 

1,J) *) 

PRINT 103, ((C (I, J) ,1 = 1,4) , J=1,M11) 

103 FORMAT(5X,F10.5,5X,F10.5,5X,F10.5,5X,F10.5) 
PRINT 99 

DO 10 1=1, NN 
READ 3, BET, DEL 

3 FORMAT ( 2F1 0 • 5) 

PRINT 9 8, BET, DEL 


C (3 y J) 


C (4 
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140 

140 

143 

144 
15 0 
153 
155 
161 
165 
165 
175 
175 
175 
177 
201 
202 

203 

204 
206 
210 
212 

214 

215 

216 
221 
223 

225 

226 
227 
230 
233 
235 
235 
237 
237 
237 
241 

244 

245 

246 
250 
252 
252 
254 
256 
26 0 
261 
264 
266 
27 0 
273 
275 
301 
304 

306 

307 
313 
315 

320 

321 
325 

325 

326 


98 FORMAT (/////5X,*8ETA =*F5.3,* DELTA =*F5.3) 
DO 11 L=1,MN 
AL=L 

S=i./(AL*BET)/DEL 

CALL SPLINE(GLAM,PHI,MM,C,S,G) 

F=G*S 

IF (AL-2.) 81,82,83 - 

81 A(i)=(l.+BET)*OEL*F 
GO TO 11 

82 A{2l = <(2.«-BET>*0EL*F-A{iJJM3. + BET) 

GO TO -11 — - - 

83 CONTINUE 
TOP=l. 

L1=L- 1 

AL1=L1 

DO 12 J=1,L1 
A J= J 

TOP=AJ*TOP 

12 CONTINUE 
L2=2*L-1 
BOT=l. 

DO 13 J=L,L2 
- -AJ=J 

EOT=(AJ+£ET)*BOT 

13 CONTINUE 
MUL=EOT/TOP 
SUM=0.0 

DO 14 N=1,L1 
AN=N 


85 

86 
87 


IF ( AN-2 •) 85,86,87 
T 0 C— 1 • 

GO TO 88 
TOD=ALi 
GO TO 88 
CONTINUE 
TOD=l. 

ICH=Ll-tN-2) 

DO 15 J=ICH,L1 


AJ = J 

TOG=AJ*TCD 

15 CONTINUE 
88 CONTINUE 

BOO=l. 

JA=L1+N 

DO 16 J=L,JA 

AJ=J 

8OD=EO0MAJ + B£T> 

16 CONTINUE 
CO=TOD/ BOD 
SUM=SUM+CO*A(N) 

14 CONTINUE 

A(L)=MUL* (OEL*F-SUM) 

11 CONTINUE 

CALL JACSERtOEL, A,EET) 
CALL NAMPLT 

CALL QIKSET(6.0, 0.0, 0. 
CALL QIKSAX ( 3,31 
CALL QIKFLT CTT,BK,101) 
CALL ENOFLT 
10 CONTINUE 
999 CONTINUE 
RETURN 
END 


0 , 6 * 0 , 0 * 0 , 0 * 0 ) 


6 

6 

6 

6 

7 

10 

11 

12 

14 

24 


SUBROUTINE JACSER (0, C,B) 
DIMENSION CC50) ,SF(50) ,P(50I 
DIMENSION 3K(101) ,TT (101) 
CQMMON/2/TI ,TF» DT,MN , BK,TT 
TT(i)=0.0 
BK ( 1) =0 . 0 
LM= 1 
T=TI 

12 T=T+OT 

X=2.*EXP (-D*T)-1. 

CALL JACOBI(MN,X»B,P) 
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26 SF ( 1) =G (1) *P (1) 

32 DO 10 L=2 ,MN 

33 L 1=L- 1 

35 AL=L 

36 SF(L>=SF(L1> +CCL)*PCL) 

43 10 CONTINUE 

45 PRINT 97,T,X • - - 

55 97 FORMATC/////5X,* T =*F6.3,* X =*F10.51 

55 PRINT 96 

61 96 FORMATC///5X,* I CCI) +,5X,* N F(T) *1 

61 DO 11 1=1,6 

65 PRINT 95, 1, CCD ,I,SFCD . 

105 95 FORMATC5X,I2,F10.2,5X,I2,F10.5) 

105 11 CONTINUE 

111 LM=LM+1 

113 BKCLM)=SFC5> 

115 TT(LM)=T 

117 IFIT.LE.TF) GO TO 12 

121 RETURN 

122 END 


7 

7 

10 

12 

14 

14 

16 

21 

22 

23 


31 

33 

34 
36 
40 

42 

43 
46 
51 
56 
64 
71 

102 

103 


C 

C 


SUBROUTINE JACOBI CN , X, E,FB> 

THIS PROGRAM CALCULATES JACOBI POLYNOMIALS OF ORDER 
K-l WITH ARG X AND PARAMETER B GT -1 
DIMENSION P B CN) 

AN=N 

IF CAN-2.) 1, 2,3 

1 PB(1)=1. 

RETURN 

2 pec i) «i. 

PB(2)=X-BMl.-X)/2. 

RETURN 

3 BSQ=B*B 
80NE=B+1. 

PB( 1) =1. 

P3C2)=X-S*Cl.-X)/2. 

DO 4 K= 3 , N 
AK=K 

AK1=AK-1. 

AK2=AK-2. 

K1=K-1 

K2=K-2 

C01=C(2.*AK1)+B)*X 

COi=( C2.*AK2)+S)*C01 

C01 = C C2.*AK2> + 90NE)* CCOi-BSG) 

C02=2.*AK2+CAK2 4-3) *( C2.*AK1) + B> 

C0=2.*AK1* (AK1+B>*<( 2.*AK2> +B) 

4 P8<K)=(CC1*PB(K1)-C02*PE (K2) >/CO 
RETURN 

END 


11 

11 


13 

14 

10 

15 

15 

11 

21 

23 

12 

23 

23 

25 

13 

26 

26 

2 

32 

34 

14 

35 

35 

3 

41 

43 

15 

43 

43 

4 


SUBROUTINE SPLINE t X , Y ,M , C , XINT , YINT > 

DIMENSION XC50) ,Y(50) ,C(4,50) 

IF C XINT-X ( 1 ) ) 1,10,11 
YI NT=Y ( 1) 

RETURN 

CONTINUE 

IFCXCM) -XI NT )1, 12, 13 
YI NT=Y( M) 

RETURN 

CONTINUE 

K=M/2 

N=M 

CONTINUE 

IF ( X C K) -XINT) 3,14,5 
YINT=YCK) 

RETURN 

CONTINUE 

IFCXINT-XCK+11) 4,15,7 
YINT=Y ( K+l) 

RETURN 

CONTINUE 

YINT=CXCK + 1)-XINT> *CCC1,K»*(XC K+l) -XINT) **2 + CC3,K>> 
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54 


65 


65 

5 

65 


70 

6 

72 


72 

—16 

74 


75 

17 

77 


— 100 

— 

100 

7 

102 


103 

8 

-103 

— 

106 

18 

106 


111 

19 

113 - 

- - - 

114 


115 

1 

121 

101 

121 


123 



7 


7 


7 


11 


12 


15 


20 

2 

26 


27 

3 

34 


37 


41 


44 


50 


51 


53 


54 


61 


65 


70 

4 

74 


76 


101 


105 


112 


114 


116 


117 


120 

6 

127 


133 


135 


140 


143 


146 


154 

7 

165 


165 



YINT=YINT*(XINT-X(K) )*(C(2»K)*(XINT-X(K))**2+C(4, K) ) 

RETURN 

CONTINUE 

IF(X(K-1)-XINT)6,16,17 
K= K-l 
GO TO 4 - 

YINT=Y (K-i) --- - 

RETURN 

N=K 

K=K/2 

-GO -TO 2 - — 

LL = K 

K= (N + K) /2 
CONTINUE 

IF(X(K)-XINT)3,14,18 - - — 

CONTINUE 

IF (X (K-l) -XINT) 6,16,19 
N=K 

K=(LL+K)/2 
GO TO 8 
PRINT 101 

FORMATt* OUT OF RANGE FOR INTERPOLATION *) 

STOP 

END 


SUBROUTINE SPLICE(X,Y,H,C) 

DIMENSION X (50) * Y ( 50 ) ,0(5 0) ,P(5Q> ,E(5Q> ,C(4,50) 
DIMENSION A (50, 3) , S( 5u) ,Z (5 0) 

MM=M-1 — 

00 2 K=1,MM 
DCK)=X(K+1) -X(K) 

P(K)=D(K)/6. 

E(K)=(Y (K+1)-Y(K))/0(K) 

DO 3 K=2,MM 
B(K)=E(K)-E(K-1) 

A(l,2) = -1.-D(l)/D(2) 

A (1,31=0 (1) /D (2) 

A(2,3)=P(2)-P(1)*A(1»3) 

A(2»2)-2«MP(1> + 3 ( 2) )-P(ll*A(i,2) 
A(2,3)=A(2,3)/A(2,2) 

8( 2)=3( 2)/A(2,2) 

DO 4 K=3,MM 

A(K,2)=2.*(P(K-1)+P(K))-F(K-1)*A(K-1,3) 
B(K)=B(K)-P (K-l)^B(K-l) 

A(K,3)=P(K)/A(K,2) 

B(K)-B(K)/A(K,2) 

Q=D(M-2)/D (M-l) 

A(M,l)=l.+Q+A(M-2,3) 

A(M,2)=-Q-A(M,1)*A(M-1,3) 

B(MJ=B(M-2)-A(M*i)*B(M-l) 

Z(M)=B(M)/A(H,2) 

MN= M-2 
DO 6 1=1, MN 

Z(K)=B(K)-A( K, 3) ♦ Z (K +1) 

Z( 1) =-A (1,2)*Z(2)-A(1,3)*Z(3) 

DO 7 K=1,MM 
Q=1 •/ (6 • *D ( K) ) 

C(1,K)=Z(K ) *Q 

C(2,K)=Z(K+i)*Q 

C(3,K) = Y(K)/D (K) -Z(K)*P(K) 

C(4,K)=Y(K+1)/0(K)-Z(K+1)*P(K) 

RETURN 

END 


40 


TofalomZ 'Unpa.ct 


PROGRAM BETA (INPUT, OUTPUT, PUNCH, PLOT, TAPE 99=PLOT) 
3 REAL N0N(4) , F (4 ,4 , 1) ,G(4,4 ) , Q( 4) ,PT(4) 

3 REAL B(4 ) ,C(4) 

3 REAL LPC50) , OTA (50) 

3 EQUIVALENCE (NON, 3) 

3 COMMON K1,K2.K3,K4 

3 C OMMON/A UX/H,P,PK1,PK2,BMU,X,Y 

3 LP(i)=0-0 

4 DTA(i)=0.0 

5. READ 2»Ki,K2*K3,K4 

20 2 FORMAT (12) 

* Ki a OROER OF SYSTEM OF EQUATIONS 

* K2 = NO. OF DISTINCT KERNELS 

* K3 = NO. OF DATA POINTS 

* K4 = NO. OF DATA SETS TO BE EVALUATED 

* SET UP DATA POINTS 


20 


A K=K3 

22 


DO 5 N=i,K3 

23 


A N=N 

24 

5 

PT(N)=AN/AK 


* SET 

1 UP INTEGRATION MATRIX 

31 


M=K3-2 

33 


N=K3-1 

34 


A =K3 

35 


A=1./(3.*A) 

37 


DO 10 K=2*M, 2 

41 

10 

0 (K)=2.*A 

46 


00 15 K= 1 , N , 2 

47 

15 

D(K)=4.*A 

54 


D (K3)=A 


* CALCULATE NONHOMOGEN EOUS TERMS 

56 


RHS=1.0 

57 


DO 22 1=1, K2 

€1 


PRINT 9 

64 

9 

FORMAT (1H1) 

64 


READ 61, 8MU 

72 

61 

FORMAT (F10'« 5) 

72 


DO 999 11=1, K4 

74 


D C 35 N=1,K3 

75 

35 

NON(N)=RHS*PT(N)*PT(N> 


* CALCULATE KERNEL MATRICES 

102 


CALL CONST ( I ) 

104 


DO 20 N= 1 , K 3 

106 


DO 20 M=1,K3 

107 


IF(M-N)25 ,30,30 

112 

25 

F (M,N,I)=F(N,M,I) 

121 


G C TO 20 

121 

30 

F (M,N,I) =FU(I,PT(H) ,PT(N) ) 

132 

20 

CONTINUE 

137 


CALL CHA NGE (F,G,D,I) 

142 


CALL LINEQ ( G, 3, C, K3) 

145 


DO 40 L=1,K3 

147 


PRINT 6, PT ( L) ,NON (L) 

156 

6 

FORM AT (5 X,F8. 4, FI 5. 6) 

156 

40 

CONTINUE 

161 


LP(II«-1) =N0N(K3) 

163 


DTA ( II+l ) =P 

165 

999 • 

CONTINUE 

1 67 


PUNCH 66, (DTA(IX) ,LP(IX) ,IX 

2 (3 

66 

F ORMAT (2F10 • 5 ) 

2 03 


CALL LAP IN V (DTA,LP) 

205 

22 

CONTINUE 

210 


END 


6 


FUNCTION SIMP (I , A , B) 

C OMMON/A UX/H, P, P Kl ,PK2,BMU , 

6 


DEL=0.25*(B-A) 

10 


IF (DEL)40,45,50 

12 

45 

S I MP = 0. 0 

13 


RETURN 

14 

50 

CONTINUE 

14 


SA=Z(I,A)+Z (1,3) 

16 


SB=Z(I,A+2.*OEL) 

T 5 


SC=Z (I,A+DEL) +Z(I , A+3.*DEL) 
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53 


Sl=(0EL/3.) MSA+2.*SB+4 

61 


IFCS1.EQ.0. 0) GO TO 45 

62 


K =3 

63 

35 

SB=SB+SC 

65 


D EL=0.5*DEL 

67 


SC=Z (I,A + OEL) 

. _75 


_J=K-1 

77 


DO 5 N=3 , J, 2 

100 


A N=N 

101 

5 

S C=SC«-Z(I,A-»-AN*OEL) 


113 
122 
125 
1 27 
131 
1 23 

133 

134 
136 
140 
152 
152 
162 
162 
166 
1 67 
171 
175 
202 
206 
207 


30 

25 


40 

42 

60 


70 


S2= (DEL/3.) *CSA«-2 .*SB«-4.*SC> 
DIF=ABS( CS2-SD/S1) 

ER=0.01 

IFCOIF-ER)30,25,25 

SIMP=S2 _ 

RETURN 
K=2*K 
S 1=S2 

IF (K— 2043) 35, 35,4 0 
PRINT 42 , I, A , 8 

FCRMATC5X,* INT. DOES NOT CONVERGE 
PRINT 60 , X, Y 
FORMAT (2F1Q»5) 

DO 70 J=1,1Q 
DIP=J 

D IP=DIP/10. 

W=Z ( I,DIP ) 

PRINT 60, W 
CONTINUE 
CALL EXIT 
END 


*,I3,2F9.4) 


7 

7 

7 

10 

11 

24 

30 

31 

40 

41 


10 


20 


SUBROUTINE CHANGE (F,G,D, I) 

REAL F(4,4,l) »G (4 ,4) ,0 (4) 

COMMON K1,K2,K3,K4 

DO 10 N=1,K3 

DC 10 M=i,K3 

G CM, N) =F(M,N,I) *D(N) 

CONTINUE 

00 20 N=1,K3 

G 1N,N)=G(N,N)«-1.0 

RETURN 

E ND 


7 

7 

10 

17 

20 

22 

23 

33 

34 
41 
43 
-44 

55 

€1 

65 

66 
10 
11 

1 -5 

1 10 
1 11 
1 14 
1 16 
1 11 
1 <2 
1 2 
1 16 


SUBROUTINE LINEQC A ,3,T,N) 
REAL A <N , N) , 3 (N ) , T (N> 

DO 5 1=2, N 

5 A (I » 1 )=A (1,1 ) /A (1 ,1) 

00 10 K=2,N 
M=K-1 

DO 15 1=1, N 
15 TCI) =A (I , K) 

DC 20 J=1,M 
A (J,K)=T (J) 

J1=J+1 

DC 20 I=J1,N 
T (I)=T(I)-A(I,J)*A(J,K) 

20 CONTINUE 

A (K,K)=T(K) 

IF(K.EQ.N) GO TO 10 
M=K *-l 

DO 25 I=M,N 
25 A CI,K)=T (I)/A(K»K) 

10 CONTINUE 
* BACK SUBSTITUTE 
DO 30 1=1, N 
T (I ) =B(I ) 

M=I + 1 

IF(M.GT.N) GO TO 30 
DO 30 J=M,N 
B CJ)=8(J)-A(J,I)*T(I) 

30 CONTINUE 

DO 35 1=1, N 
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. 27 
l tl 
.46 
.50 
L 51 
.52 
L £4 
.62 
L 67 
.67 


K=N+ 1 -I 

B <K)=T(K)/A (K» K) 
K 1 =K -1 

IF(Kl.EO.O) GO TO 35 
DO 35 Jl=l,Ki 
J=K-Ji 

T (J)=T(J)-A(J,K)*B CK) 
35 CONTINUE 
RETURN 
END 


6 

6 

7 

10 

11 

12 

13 

;o 

21 

23 

25 

22 

23 

‘6 

27 

U 

17 

t -7 


10 

5 

20 


15 


FUNCTION FU(I>A, 3 ) 
COMMON/AUX/H,P,PKi,PK 2 ,BMU,X,Y 
X=fl 
Y =B 

IF(A* 3 ) 5 , 10,5 
F U= fl • 0 
RETURN 

SUM-SIMP ( 1 , 0 . 0 , 5 , 0 ) 

ER-OaOl 
DEL = 5.0 . 

UP=DEL+ 5.0 

ADDL=SIMP(I,D£L,UP> 

DEL =UP 

TEST=ABS (ADOL/SUM) 

SUM=SUM+ ADOL 

IF C TEST-ER) 15 , 20 , 20 

FU=SQRT(X*Y)*SUM 

RETURN 

END 


3 

3 

5 

6 
15 
24 
31 
31 

33 

34 

36 

37 

41 

42 
57 


57 

60 


SUBROUTINE CONST (II 
COMMON/AUX/H,P,PK 1 ,PK 2 , 8 MU,X,Y 
PR 1 = 0 .Z 9 
PR 2 = 0.29 

PK 1 =SQRT ( (i.- 2 .*PRl)/( 2 .Ml.-PRl>n 
PK 2 = SQRT ( (l.- 2 .»PR 2 )/( 2 .*(l.-PR 2 ))l 
READ 1 ,P 

1 F ORMAT (FIG ■ 5 ) 

HF= 5.0 

HH= 0 .2 
H H= 0 .5 
H H=i . 0 
H F= 2.0 
H-l. /HH 

PRINT 2 *BMU*PR 1 ,PRZ,HH,P 

2 FORM AT<///// 5 X, * MU 2 /MU 1 =*F 6 . 2 ,* NU 1 =*F 4 .Z,* NU 2 =*F 4 . 2 /// 5 X,* A 
1 /H =*F 4 . 2 ,* C 21 /PA =*F 4 . 2 ) 

RETURN 

END 


5 

5 

26 

30 

32 

33 

34 
34 
36 
45 
54 
60 
62 
72 
.02 
.05 
.16 

17 


FUNCTION Z(I*S> 

COMMON/AUX/H,P,PKl,PK 2 ,BMU,X,Y 
BESJT <A)=SQRT(Z.*A/PIIMSIN(A) /A/A-COS (A > / A. I 
PI= 3. 1415926 
IF(S- 0.0 ) 5 , 5 , 10 
5 Z = 0 . 0 
R ETURN 
10 C CNTINUE 

PP=P*P 

GB=SQRT(S*S+ 1 ./PP) 

GD=SQRT ( S*S+ 1 •/ BMU/PPJ 

AA= 1 .- 3 MU*GD/G 3 

AB=i.+BKU*GD/G 8 

AC= 1 .-AA/AB*EXP(- 2 .*GB*H) 

AD= 1 .+AA/AB*EXP(- 2 .*GB*H) 

F=G 3 *AC/A 0 

Z=(F-S)* 3 ESJT (S*X ) *BESJT (S*Y) 

RETURN 

END 


43 - 



c 

c 

c 


SUBROUTINE LAPINV (GLAM, PHI) 

THIS PROGRAM EVALUATES THE COEFFICIENTS FOR SERIES 
OF JACOBI POLYNOMIALS WHICH REPRESENTS A LAPLACE 
INTEGRAL 


5 

-3 


5 

5 

5 


16 

1 

16 

30 

2 

30 

34 

.99 

34 

40 

44 

101 

44 

65 

102 

65 

67 

73 

303 

73 

112 

103 

112 
1 16 
1 cl 
1 30 

3 

130 

140 

98 

140 

143 

144 
150 
153 
1 55 
1 Cl 

81 

1 C5 
1 '5 

82 

1 "5 
1 *5 

83 

1 75 

1 77 

2 fl 
2 IZ 
2 C3 
2 T4 
2 C6 

12 

2 10 
2 12 
2 14 
2 15 
2 16 
2 21 

13 

2 -3 
2 25 
2 -6 
2 c7 
2 20 
2 23 

85 

2 25 
2 5 

36 

2 27 
2 27 

87 

2 *7 
2 *1 
2 *4 
2 5 
2 *6 
2 CO 

15 

2 '2 

83 

2 C2 
2 24 
2 r 6 
2 CQ 



A (50)., GLAM (50) ,PHI (5 0) ,C (4,50) 
BK(iOl) * TT ( 101 ) 

• MN,BK, TT 


PHI 


C (2 , J) 


*> 


INVERSION 
REAL MUL 
DIMENSION 
DIMENSION 
C0MM0N/2/TI*TF, OT 
REAO 1,NN,MN,MM 

FORMAT (312) ...... 

READ 2,TI,TF,DT 
FORMAT (3 F1Q .5) 

PRINT 99 

FORMAT (1H1) _ 

CALL SPL ICE (GLAM, PHI, MM, C) 

PRINT 101 

FORMAT (/////5X»* GLAM 
PRINT 102, (GLAM (I ) , PHI (I) ,1=1, MM ) 

FORMAT (5X,F10*5,5X,F10*5) 

M 11= MM-1 
PRINT 300 

FORHAT (////5X Cfl,J) 

.♦J) *) 

PRINT 103, ((C(I,J) ,1=1,4), J=1*M11) 

F CRM AT (5X,F10 .5 ,5 X ,F10,5 ,5X , F10 . 5, 5X ,F10 . 5 ) 
PRINT 99 
DO 10 1=1 ♦ NN 
READ 3, BET ,OEL 
F ORMAT (2 F10 • 5 ) 

PRINT 98, BET, DEL 

FORMAT (/////5X**B ETA =*F5.3,* DELTA =*F5.3) 
DO 11 L=l, MN 
A L=L 

S=1 * / (AL + 3ET) /DEL 

CALL SPLINE (GLAM* PHI, MM,C, S,G) 

F=G*S 

I F(AL-2a ) 81, 82, 83 
A (1 ) = (I* ♦BET) *DEL*F 
GO TO 11 

A (2)=((2.+BET)*0EL*FrA(l))*(3.+B ET) 

GO TO 11 
CONTINUE 
T OP = i a 
L1=L-1 
A Li=Ll 

DO 12 J=1,L1 
A *=J 

TOP=AJ*TOP 

CONTINUE 

L2=2*L-1 

BOT=l. 

DO 13 J=L,L2 
A J=J 

BOT= (AJ+3ET) *30T 

CONTINUE 

MUL=BOT/TOP 

SL’M=0.0 

DO 14 N=l,Li 

A K=N 

IF ( AN— 2# )85,86, 87 
T 0 D= 1 # 

GO TO 88 
T0D=AL1 
GO TO 83 
CONTINUE 
TOD=l. 

I CH=L1- ( N-2 ) 

DO 15 J=ICH,Li 
A J= J 

TOD=AJ*TOD 
CONTINUE 
CONTINUE 
B 00=1. 

JA=L1*N 
DO 16 J=L» JA. 

A J=J 


C(3, J) 


C (4 
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2 ei 

800= BOD* (AJ + BET ) 



2 64 

16 CONTINUE 



266 

C 0=T OD/BOD 



270 

SUM=SUM*CO*A(N) 



273 

14 CONTINUE 



275 

A (L)=MUL*(OEL*F-SUM) 



301 

11 CONTINUE — 


. . — 

304 

CALL JACSER (DEL , A , BET ) 



306 

CALL NAMPLT 



307 

CALL QIKSET (6.0 ,0 • C,0.0,6.0,0.0,0.0) 



-_.313 

... CALL-QIKSAX (3,3) 

. . 

... . . 

315 

CALL QIKPLT (TT,BK,1Q1) 



3 20 

CALL ENDPLT 



3 21 

10 CONTINUE 



3 25 

999 CONTINUE 




3 25 

RETURN 



3 26 

END 




SUBROUTINE JACSER (0,C,B) 


6 


DIMENSION C(5Q) ,SF (50),P (50) 


6 


DIMENSION BK(101) ,TT(10i) 


6 


C0MM0N/2/TI, TF, OT , MN, BK,TT 


6 


TT(1)=0.0 


7 


BK ( 1 ) =0 • 0 


10 


LM=1 


31 


T =T I 


12 

12 

T=T*OT 


14 


X=2.*EXP(-D*T)-1. 


24 


CALL JACOBI (MN* X, 3»P) 


26 


SF(l)=C(l)*Pll> 


32 


DO 10 L= 2 , MN 


33 


L 1-L-l 


35 


A L=L 


36 


SF(L)=SF(L1) ♦CIL) *P(L) 


43 

10 

CONTINUE 


45 


PRINT 97,T,X 


55 

97 

F0RMAT(/////5X,* T =*F6.3,* .X =*F10.5) 


55 


PRINT 96 


61 

96 

FORMAT (///5 X,* I C(I) *,5X,* N F (T ) 

*) 

61 


DO 11 1=1,6 


65 


P RINT 95 ,I»C ( I) ,1 ,SF(I) 


1 05 

95 

FORMAT (5X*I2»FIQ, 2 ,5X,I2,F10 .5 ) 


105 

11 

CONTINUE 


111 


LM=LM+1 


113 


BK(LM)=SF (5) 


115 


TT(LM)=T 


1 17 


IF(T.LE.TF) GO TO 12 


1 21 


RETURN 


1 22 


END 



SUBROUTINE JACO 31 C N, X , 3, P3) 



C 


THIS PROGRAM CALCULATES JACOBI 

POLYNOMIALS OF OROER 


C 


K-l WITH ARG X AND PARAMETER 8 

GT -1 

7 



D 1MENSIO N PB ( N) 


7 



A )=N 


10 



I F(AN-2. )1,2,3 


12 


1 

PB(1)=1. 


14 



RETURN 


34 


2 

P B ( 1 ) =i« 


16 



PB(2)=X-B*(l.-X)/2. 


21 



RETURN 


2 


3 

BSQ=B*B 


23 



30NE=B*1. 


25 



PB(1)=1. 


26 



PB(2)=X-B*(1.-X)/ 2. 


:± 



DO 4 K=3 ,N 


3 



A K=K 


•4 



AK1=AK-1. 


26 



A K2 = AK— 2 • 


iO 



Ki=K-i 


-2 



K2=K-2 


<3 



C01= ((2.*AKi) +3 )* X 


i 6 



C01= ( (2. *AK2) +B> *C01 
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51 

56 

64 

71 

02 

03 


C01= ( (2. *AK2) +BON E )* (COi-BSQ) 
C02=2.*AK2* (AK2 + B )*( C2.*AK1H-B) 
C0=2.*AK1*(AK1+B> * ( (2.*AK2) + 3) 
PB(K)=(C01*PB(K1) -C02*PB(K2) )/CO 
RETURN 
END 


SUBROUTINE SPLINE t X, Y,M, C, XINT,YINT> 

( 50 ) ,C (4, 50 ) 

*11 


13 


11 


DIMENSION X ( 5 0) ,Y 

11 


IFCXINT-X(l) 11,10 

13 

10 

YINT=Y(1) 

14 


RETURN 

.15 

11 

CONTINUE 

15 


IF (X (Ml — XINT1 1, 12 

;i 

12 

YINT=Y(M> 

-3 


RETURN 

',3 

13 

CONTINUE 

23 


K=M/2 

;5 


N=M 

76 

2 

C CNTINUE 

lb 


IF ( X (K) — XINT) 3,14 

72 

14 

YINT=Y(K) 

4 


RETURN 

75 

3 

CONTINUE 

75 


I F { X I NT- XlK+1114, 


15 

YINT=Y(K+1J 

IZ 


RETURN 

LZ 

4 

CONTINUE 

i 3 


YINT= (X ( K+I) “XINT 

54 


Y INT=YINT* ( X INT-X 

65 


RETURN 

65 

5 

CONTINUE 

65 


IF(X (K-l)-XINT) 6, 

70 

6 

K=K-i 

72 


GO TO 4 

72 

16 

YINT=Y (K—11 

74 


RETURN 

75 

17 

N=K 

77 


K=K/2 

100 


GC TO 2 

100 

7 

L L=K 

102 


K=(N+K)/2 

103 

8 

CONTINUE 

103 


I F(X (K)-XINT) 3,14 

106 

18 

CONTINUE 

106 


I F ( X ( K“1 1 “X I NT) 6, 

111 

19 

N=K 

113 


K = (LL+K) / 2 

114 


GC TO 8 

1 15 

1 

PRINT 101 

1 21 

101 

FORMATS OUT OF R 

1 21 


STOP 

123 


END 


15,7 


*{C(1,K)MX{K*1)-XINT) **2 + Ct3 ' 


13 


FOR INTERPOLATION *1 


7 

7 

7 

11 

12 

15 

20 

26 

27 

34 

37 

41 

44 

50 

51 

53 

54 
61 


SUBROUTINE SPLICE < X,Y,M, Cl 

DIMENSION XC50),Y (50 1 ,0 (50 1 , P( 50 1 , E ( 50 1 , C C 4, 50 1 
DIMENSION A(50*3) ,B(50),Z(50) 

MM=M-1 

DO 2 K=1 • MM 

D(K)=X(K+1)-X{K) 

P (K)=D(K)/6. 

E <K) = (Y(K*1)-Y(K) 1 /O(K) 

DO 3 K=2 * MM 
B (K)=E(K)“E(K-1) 

A <l,2)=-i.-D(i) /D C 2) 

A <1 » 3)=D (1)/D (2) 

A (2, 3)=P (2) -P(ll* A (1,3) 

A C2,2)=2.*(P(1) *P (2) >-P(l)*A(l,2> 

A (2, 3)=A (2, 3) /A(2 ,2) 

B (21=3(2) /A (2,2) 

DO 4 K=3 , MM 

A (K,2)=2.MP(K-1) +P(K))-P(K-1)*A (K-1,3) 

B (K)=B(K)-P(K-i>*B(K-l) 
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65 
70 
74 
76 
101 
105 
112 
114 
1 16 
117 
J.20 
127 
1 23 
135 
140 
143 
146 
154 
165 
165 


A (K,3) = P(K) /A (K *2 ) 

4 3(K)=B(K)/A (K,2> 

Q=0 (M-2) /OlM-1) 
A(M,i)=l.+Q+A(M-2,3> 

A (M,2)=-Q-A (H,l)*A(H-i,3) 
3 (M>=B<M-2>-AtM,l)*B(M-l) 
ZCM)=BCM)/AIM,2> 

MN=H-2 
DO 6 1=1, MN 
K=H-I 


_6 Z(K)=B(K)rA(K,3)*Z(K+i) 

Z(1)=-A(1,2)*Z(2) -A Cl, 3) *Z ( 3 ) 
00 7 K=1 ,MM 
Q=l./(6. *D<K) > 

C (1,K)=Z (K) *Q 
C C2,K)=Z CK+1) *Q 
C (3,K)=Y (K)/DCK)-Z(K)*P(K) 

7 C (4,K)=Y (K*l>/0(K)-Z(Kfl>*P(K> 
RETURN 
ENO 
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